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We develop a new analytic approach for the study of lattice heteropolymers, and apply it to 
copolymers with correlated Markovian sequences. According to our analysis, heteropolymers 
present three different dense phases depending upon the temperature, the nature of the 
monomer interactions, and the sequence correlations: (i) a liquid phase, (ii) a "soft glass" 
phase, and (m) a "frozen glass" phase. The presence of the new intermediate "soft glass" 
phase is predicted for instance in the case of polyampholytes with sequences that favor the 
alternation of monomers. 

Our approach is based on the cavity method, a refined Bethe Peierls approximation 
adapted to frustrated systems. It amounts to a mean field treatment in which the nearest 
neighbor correlations, which are crucial in the dense phases of heteropolymers, are handled 
exactly. This approach is powerful and versatile, it can be improved systematically and 
generalized to other polymeric systems. 

PACS numbers: 81.05.Lg, 64.70.Pf, 36.20.Ey 



I. INTRODUCTION 
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In the last 20 years much effort has been devoted to the theoretical study of heteropolymers 



_ 58 . One of the main motivations was to understand the statistical physics of protein folding 
J 3> 12) Q Ifi^ ]- Despite the insight that has been accumulated, the goal remains distant. 
On the one hand, most analytical studies have been limited to random bond models j^ . 6f] | (in 
which the interaction energies of all the couples of monomers alon g th e chain are independent 
random variables), or to uncorrelated random copolymer sequences jigl. However, there are 
many indications that sequence correlations induced by natural selection play an important role 
for the folding and stability of proteins. On the other hand, in this difficult problem, analytic 
computations have to resort to some approximations which are not easy to control. It is thus 
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important to have a variety of different techniques at hand in order to crosscheck the predictions. 

In this paper we develop a new tool for the analytical study of heteropolymers, based on the 
cavity method as used in various frustrated systems (a short account of our results has appeared 
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43j). We use this method to investigate the phase diagram of copolymers with Markovian 



sequences. Within our approach we find copolymers to exist in three distinct dense phases (apart 
from the diluted coil phase at high temperature) depending upon the structure of the interaction 
energy matrix, the sequence correlations and the temperature: (i) The liquid globule phase in 
which distinct monomers are essentially uncorrelated and can freely rearrange within the globule 
(apart from obvious constraints on monomers that are close along the chain); (ii) the "frozen 
glass" phase in which the polymer is stuck in one out of a few well-separated low-energy confor- 
mations; (Hi) a "soft glass" phase with broken ergodicity (in the thermodynamic limit) in which 
the thermodynamically relevant conformations form a continuum in configuration space. This last 
phase has never been predicted in an analytical computation (although such a possibility has been 



and a very similar phase seems to be present in 



envisioned in phenomenological models 
the numerical results of on the dynamics of heteropolymers.). Albeit frustrated, it has a much 
larger entropy, and appears already at a smaller density than the usual "frozen glass" phase. 

Some of the most successful tools used so far in the study of random heteropolymers are mean 
field approaches based on the replica method 0,0, Q. Crucial to these calculations was the 
identification of some relevant order parameter, and the proposition of a suitable Ansatz describing 
the phase transition in a coupled space of real space coordinates and replica indices. This type 
of approach is potentially very powerful, but it becomes quite complex for heteropolymers. On 
the one hand, it requires a physical intuition for identifying the relevant degrees of freedom and of 
their behavior. On the other hand, an Ansatz tailored to describe a certain type of physics may 
hide other, unexpected features. 

Our cavity method consists in a refined version of the Bethe Peierls approximation. While 
this also represents a kind of mean-field approximation, it differs fundamentally from the previous 
ones. Applying the Bethe-Peierls approximation to lattice heteropolymers allows to describe self- 
consistently the frustration on a local microscopic level. This approach can be thought of as 



331 ] . Its general 



the first step in the series of cluster variational (or Kikuchi) approximations 
philosophy consists in keeping track of local correlations inside some small region exactly, while 
treating the external degrees of freedom as an environment whose statistical properties have to be 
determined self-consistently. In the Bethe approximation, the only correlations which are treated 
exactly are the ones between neighboring sites on the lattice. This is an improvement with respect 
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to the naive mean field that treats distinct sites as statistically independent. Moreover, it is the 
first of such approximations to be meaningful fo r poly mers, since the backbone structure induces 
strong correlations between neighbors Q, Q, il U, E , 42 ] • Another potential advantage of the cavity 
method is that it can be used for one given polymer, without the need to average over an ensemble 
of sequences as in the replica method. While in the present work we focus on ensemble-averaged 
properties, one should keep in mind this possibility which could lead to interesting algorithmic 
developments in the future. Finally, the refined Bethe-Peierls approximation is supposed to be 
exact on locally tree-like structures (e.g., on random graphs). This is an important feature: It 
allows one to set up the mean-field analysis in a mathematically well-defined way, and its predictions 
can be checked against numerical simulations on those random "mean-field" lattices for which the 
theory is expected to be exact. 

Within our cavity method, any heteropolymer is found to undergo a glass transition at large 
enough densities. Two main schemes of glass transitions can occur, depending on the details of 
the sequence, each of them being associated with one of the types of glasses mentioned above. 

The transition to the frozen glass phase is a discontinuous transition, which is called random 
first order, or one step replica symmetry breaking (1RSB) transition in the replica language. It 
corresponds to the type of transition which has been found in many previous studies, of which the 
Random Energy Model (REM) ^| is the simplest archetype. 

The transition to the soft glass phase is a continuous one, corresponding to full replica symmetry 
breaking (FRSB). This is more in line with recent scenarios proposing a freezing that proceeds 
gradually from small scales to larger and larger structures In a series of papers exploiting 

a Gaussian variational technique to deal with the dynamics of heteropolymers, copolymers in 



particular, a much richer phase diagram was proposed, where the ultimate REM-like folding to a 

3, El 



unique ground state is preceded by a less structured but still frustrated glassy phase 
As for the glass transition, the random copolymer was proposed to be in the same universality 
class as the Ising spin glass ^f|, which would imply a continuous transition with a full breaking of 
the replica symmetry. 

Beside providing an alternative and well controlled analytical approach, our cavity analysis adds 
to the above pictures in that it highlights the dependence of the scenario to be expected on the 
correlations of the monomer sequences. 

In order to keep the computations more transparent we avoid here the use of replicas (although 
it would be possible to write all of the ensemble-averaged cavity equations using replicas), but we 
keep to the traditional replica vocabulary of 1RSB and FRSB to denote the two types of transitions. 
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We will apply here the general method to treat Markov-correlated sequences. However, a much 
wider range of possible applications of this technique is open. 

The paper is organized as follows: In Section [n] we define the lattice model and review the 
treatment of polymers in the grand canonical ensemble. We then introduce the basic ideas of the 
Bethe approximation and discuss the G-collapse from the random coil to the liquid globule phase. 
Section ITTT1 discusses the shortcomings of the liquid solution and generalizes the method to the case 
where many pure states exist (as typically in a glassy phase). In particular, we propose a set of 
local order parameters that allow to distinguish both theoretically and experimentally between two 
different types of glass transitions. In Section HVl we describe some basic tools for analyzing the 
glass transition. We present a local stability criterion for the liquid phase and the 1RSB cavity 
equations which are used to describe the glassy phase. This formalism is illustrated in Section IVlbv 
considering the exemplary cases of alternating sequences with attractive or repulsive interactions 
of like monomers. 

It turns out that the two types of interactions imply very different phase transitions: either 
a continuously emerging "soft" glass phase or the "standard" discontinuous freezing transition. 
These two scenarios are found in the study of Markovian chains in Sec. IVII The properties of the 
strongly frozen phase is analyzed in Section fVIII bv focusing on maximally compact conformations. 
We conclude with a summary of our results and a discussion of their relevance for protein folding. 
Several technical developments are included in the seven appendices. 



II. THE CAVITY APPROACH TO HETEROPOLYMERS 



In this Section we describe the type of heteropolymer models which we shall study. We derive 
their phase diagram under the assumption that the polymer is "liquid" meaning that any sta- 
tistically relevant conformation is dynamically accessible to the molecule. In replica jargon this 
corresponds to assuming replica symmetry. The next sections will render more precise the regions 
of the phase diagram where this liquid phase is stable and corresponds to the physically relevant 
state. 



A. The lattice polymer model 

Our starting point is the standard model of lattice polymers [ill. I^. which we generalize for 
polymers living on a general graph Q. We denote by i, j, . . . £ V the vertices of Q (with |V| = V), 
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and by S 8 the edges of Q. Let to = {u\ . . . u>n), ^>a £ V denote a self-avoiding walk (SAW) 

of length N on Q. The position of a monomer along the chain is denoted by a, b, . . . £ {1 . . . N}, 
and we assume an interaction matrix e a b to be assigned. The corresponding energy reads: 

H N (uf) = ^ e ab, (1) 

(a,b)\(uj a ,uj b )G£ 

where the sum runs over couples of non-consecutive monomers which are nearest neighbors on the 
lattice. 

The choice of the matrix e a b is crucial. The standard homopolymer model is recovered by 
setting e a b = eo- A popular model in heteropolymer studies is the random bond model which 
assumes the e a b to be independent identically distributed (i.i.d) quenched random variables. In 
this work we study the more realistic case where the interaction energies are determined by the 
underlying monomer sequence. The sequence will be given by {a±, . . . ,ctn}, with a a £ A being 
the type of the monomer at position a in the sequence. The interaction energy of two monomers 
is assumed to depend only upon the monomer type: e a b = E aa0 - b . In particular, we shall focus 
on copolymers (although the approach is general) where there are only two types of monomers: 
A = {A,B}. Interaction matrices E a%a i of particular interest are: 

• The HP model. A and B monomers represent (respectively) hydrophobic and polar 
aminoacids, and the interaction matrix is chosen accor ding ly e.g., Eaa = — 1> Eab = 
Ebb = 0. This is a popular toy model for protein folding |15| |. 

• The polyampholyte. A and B are supposed to carry screened charges which suggests Eaa = 
Ebb = +1 and Eab = Eba = — 1- Sometime we shall refer to this interaction matrix as the 
antiferromagnetic (AF) model. 

• The symmetrized HP model. We take Eaa = Ebb = —1 and Eab = Eba = +1- This is 
the standard model for copolymers with monomers that have a tendency to segregate js^ . 
We shall refer to it as the ferromagnetic (F) model. 

As for the graph Q we shall consider two particular cases: (i) A T/-sites portion of the d- 
dimensional cubic lattice, (ii) A l/-sites Bethe lattice, i.e., a random lattice with connectivity 
(k + 1). Its interest stems from the observation that, in the thermodynamic limit, our mean-field 
calculations are exact on such a graph. 

Both for our analytical computations and for the simulations on the Bethe lattice we shall 
need to consider periodic sequences with period L: Oi = cri+L- The complete sequence is therefore 
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determined by its first period (a\ . . . Hereafter, we shall use the shorthand notation "monomer 
a" to refer to all monomers in positions a + nL with integer n. Furthermore, monomer indices 
always should be read modulo L. We expect the non-periodic case to be recovered in the L — > oo 
limit, even if this limit is taken after the limit N, V — > oo. 

The random-bond model is obtained in the \A\ = L —* oo limit by setting o\ 7^ 02 7^ • • • 7^ 
and taking the E aa i to be i.i.d. random variables. 

In order to understand the influence of the correlations in the sequence of monomers, we shall 
consider Markovian random copolymer chains in the large L limit. In these chains the probability 
of a monomer to be of a certain type depends only on the preceding monomer in the sequence. For 
the sake of simplicity we assume the two types of monomers to occur with the same frequencies. 
The statistical ensemble of the chains is then fully characterized by the probability ir £ [0, 1] of a 
monomer to be of the same type as the preceding one. 

We study the system at thermal equilibrium at a temperature T = 1/(3. We define a canonical 
free energy density as 



where the expectation value Eg is taken with respect to the graph ensemble (whenever Q is a 
random graph). The L —* 00 limit, and the expectation with respect to the sequence (oi . . .<jl) 
are (eventually) taken afterwards. 

The two free energies defined above satisfy the usual Legendre transform relation ul(/3,p) = 
Il((3,p) — pp. In order to describe free polymers (in equilibrium with the solvent) the chemical 
potential has to be adjusted to the critical value p c such that u>l(p c ) = 12]. In the grand- 
canonical picture this critical line corresponds to a phase transition between an infinitely diluted 
phase for p < p c and a dense phase with non- vanishing osmotic pressure for p > p c . If this 
phase transition is continuous, the density on the coexistence line vanishes, while it is finite if 
the transition is first order. On this coexistence line, the tricritical point where the nature of the 
transition changes is nothing but the 0-point where the collapse of the unconstrained polymer 
takes place. 

,n a hom „p„ lyme , the a, OT e d e scripta captures the e ss eet laJ „ f tne phfB e di a gram Q. 
However, in a heteropolymer, the low temperature dense phase will be strongly influenced by the 




N=pV 




and its grand-canonical counterpart 




(3) 
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FIG. 1: Possible conformations of a site (or oriented edge) on the regular Bcthc lattice. The cavity site is 
considered as the root of a branch with k leaves (here k = 5). The thick lines and filled circles represent the 
chain backbone and monomers. 

sequence heterogeneity. Due to the connectivity of the polymer chain it is in general impossible to 
find a compact folding where all interactions are favorable. The system is frustrated, and a glass 
transition will take place at sufficiently low temperature. 



B. The Bethe Peierls approximation 

As already mentioned, the Bethe approximation is asymptotically exact on locally tree-like 
graphs. Following j^jj, we define a Bethe lattice as a random lattice with fixed connectivity. Such 
a lattice is locally tree-like since the typical loop size diverges as log V with the lattice size. In 
order to handle the heteropolymer problem on a d-dimensional hypercubic lattice within the Bethe 
approximation, our approach idealizes the graph as a Bethe lattice with the same connectivity, 
k + l = 2d. 

The local tree structure of the graph can be exploited in a recursion procedure. Suppose for a 
moment that the lattice is a tree, and let us single out a single branch of the tree which is rooted 
at one 'cavity site' having only k neighbors i = 1, .., k. In the absence of 0, the branch would 
become a collection of k other branches, rooted at i = 1, .., k. This structure allows for a recursive 
computation of the probabilities of the polymer's conformations on the tree. 

We first list the possible local conformations of the cavity site in its branch (see Fig. (0): 
the site is unoccupied; (j a) or (j a): the site is occupied by the monomer a and the backbone 
continues towards the remainder of the tree, with monomer a — 1 or a + 1, respectively; (2a): the 
site is occupied by monomer a, but the polymer returns back to the leafs. (On a real tree the parts 
of the polymer on different branches are necessarily disconnected. However, on the Bethe lattice 
this is no longer the case and the polymer may be present on more than two leaves.) 
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For each local conformation a G {0, f a, J. a, 2a} of the root site 0, we denote by p^ the 
corresponding probability (as given by the Boltzmann measure). The (3L + 1) dimensional vector 
of weights p(°), with components pa\ can be expressed in terms of the corresponding k weight 
vectors pW on the neighboring sites. Note that p$ is the Boltzmann weight for the configuration 
a on i when the site is absent. We will refer to these weight vectors on root-sites as cavity fields. 

The mapping between cavity fields, p(°) = / [p*- 1 - 1 , ...,p( fc )] , can be written explicitly as: 

pf=C-'{[4\ (4) 

i=i 

p { $ = c-^j: P f a+1 ii^\ ( 5) 

i=l j^i 
i=l j^i 

p { £ = c-^Y,p ( iM^ II (7) 

where C = C[{pW}] is a normalization constant which enforces the condition YlaP^ = 1 and we 
have introduced the quantities 

# = pf + E & , = p? + E p2}c-^«.' . (8) 

a'=l a'=l 

The full lattice is built by merging + 1 branches. Therefore, once the cavity fields have been 
computed, one can express any local quantity using the neighboring cavity fields. The monomer 
density pW a t site i is a function of the k + 1 cavity fields p^') on the j = 1, k + 1 neighboring 
sites of i (recall that p^) gives the probability of a local conformation on j in the absence of i): 



P 



(0 



a=1 ii^j2 

where we have defined the normalization constant 



Pu+iPja-i lljftuh ^ d (!) , 
" «,«(p(i),..,p(*+i)) ' 



fc+l L 



^( P (i) ,-,p (fe+i) ) = n^o+^E e c-i«i n ( io ) 



The internal energy of a link (i, j) can be written in terms of the cavity fields on i and j 
(giving the probabilities of local conformations on i and j in the absence of the link 

o,6=l ^ W; (pW,pW) 
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where re^ (a, b) is the probability of having a contact between two monomers a and b along the link 
(ij) of the graph. The normalization wi(p^' ,p^') is given by 

ufV^ P W) = E iSS^ + j + & w + pSp^) .(i2) 

a,b=l o=l 

For each edge of a given graph, one can introduce a pair of cavity fields, describing 

respectively the probability of local configurations of the two points i and j in the absence of the 
edge One can write a Bethe free energy, which is a functional of all these cavity fields and 

has Eqs. ©-0 as stationarity conditions. It reads 

V(3u[{p®}) =-J2 log[^ } ] + £ log[w^] , (13) 

iev (ij)es 

where Ws and are the expressions given in and (j!2j) . respectively. Notice moreover 

that the density Q and the internal energy (jllj) can be obtained by differentiating the Bethe free 
energy with respect to the chemical potential fi and the inverse temperature f3. 

It is easy to show that the above expressions are exact if the graph Q is a tree. On a general 
lattice it holds approximately to the extent that one can neglect the correlations between the fields 
on the k + 1 neighbors of any site i, once the site i itself has been deleted. 

On a Bethe lattice, since the typical loop size diverges as log V in the large- limit, these 
k + 1 sites neighbors of i are generically distant from each other, when i is absent. Therefore 
the correlations of their fields can be beglected, if the system is in a single pure state: at low 
temperature the Gibbs measure usually has to be decomposed into pure states, within which the 
correlations between two sites decay with their distance along the graph. We thus expect the above 
cavity approximation to become asymptotically exact, insofar as cavity fields are computed within 
one pure state. 

C. The liquid solution and the 0-point 

Both on the random Bethe lattice and on the (i-dimensional cubic graph, each site has generically 
the same environment within any distance R (as long as R is kept finite in the V — > 00 limit). A 
liquid phase is therefore expected to enjoy translational invariance and will be described by a set of 
fields pa that is independent of the site. We thus look for a fixed point Pa = p* of the recursions 

©-CD- 
it turns out that the liquid solutions can be found by solving a system of \A\ +2 non- linear 

equations, |^4| being the number of monomer species in the model. This is a great complexity 
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reduction with respect to the 2>L + 1 equations @-(0). The task can be further simplified by 
using particular symmetries of the interaction matrix. This is, for instance, the case of the F- and 
AF-models defined in Sec. Ill Al which are symmetric under the interchange A <-> B. We refer to 
App. 0for a detailed discussion of how the solution is obtained. 

As shown in AppendixEIall the thermodynamic quantities depend upon the sequence {a\ . . . a£) 
only through the fractions v a of monomers of type a. As a byproduct, the L — > oo limit can 
immediately be taken. The physical meaning of this result is easily understood. In the liquid 
phase, the correlations induced by the sequence play some role just along the chain, and their net 
effect vanishes at large distance. In particular, the monomer a is surrounded by a certain fraction 
of monomers of type a' which only depends on the type of a, a a (apart from the sites occupied by 
the monomers a — 1 and a + 1, of course). 

Let us now discuss the various solutions of liquid type. 

The random coil phase is described by the trivial solution p* a = S a Q, which exists for any 
choice of the parameters. This phase has vanishing grand potential oj and density p. At high 
temperatures this is the only solution when \x is smaller than the critical chemical potential /i c 
given by exp(/5/i c ) = 1/k. At /i c a non-trivial solution emerges continuously. The latter describes 
a liquid phase under pressure (u> > for /j, > fj> c ) with a density that vanishes on approaching the 
critical line. 

The collapse of a free polymer from the random coil state to the liquid globule occurs at the 
so-called ©-point. In the grand-canonical description, it appears as the tricritical point on the line 
exp(Pp) = 1/k. Expanding around p* a = S a ^, one obtains the following relation which determines 
the ©-point temperature 



see App. ^ This result has previously been obtained within the framework of the standard cluster 
variational method j^. At temperatures below the ©-point, > (3q, the grand-canonical phase 
transition becomes first order (see Fig. I^J). The critical line ^ c (/5) is obtained by equating the 
grand potentials in the coil and globule phases, i.e., by solving oj = for the globule solution. The 
density, internal energy, and free energy are obtained by plugging the globule solution p* a into Eqs. 



In the low temperature region (3 > /3q, the dense solution can be continued to values of the 
chemical potentials smaller than the critical one // c (/3), and ceases to exist on a spinodal line. 
Likewise, the trivial dilute solution stays locally stable beyond the coexistence line up to the 



a,T£A 




(14) 
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FIG. 2: The phase diagram corresponding to the liquid (translation invariant) solution in the grand canonical 
ensemble. Above the 0-temperature, < (3q, the phase transition from the random coil phase (/i < /i c ) to 
the globule solution with finite density (/i > fj, c ) is continuous. At low temperatures, (3 > (3q , the transition 
becomes first order and is accompanied by two spinodals. The globule solution on the critical line describes 
a free polymer in coexistence with the surrounding pure solvent. The free polymer undergoes a collapse 
transition at the 0-point 



spinodal exp(/3/i) = 1/k. 

The above results compare reasonably with the outcomes of numerical simulations on d- 
dimensional lattice. For instance, the homopolymer O-point on the cubic lattice given by Tq = 1.50 

n n n n 

for d = 2 |24|, 3.716(7) for d = 3 [6J], and 5.98(6) (d = 4) |5J. Moreover the authors of Ref. |35j 
found T® = 2.25(10) on the three-dimensional diamond lattice (connectivity k + 1 = 4). These 
results should be compared with the outcome of the Bethe approximation, cf. 1)141) . which yields 
7e,Bcthc « 2.4663035 (for k = 3), 3.4760595 (k = 4), 4.4814201 (k = 5). 6.4871592 (k = 7). As for 
heteropolymers, the authors of Ref. 2^13] estimated T® 1.2 both for the F- and AF-models of 
Sec. Ill Al in d = 3. This result is compatible with Tq = l/log(2) ~ 1.442695 which comes out of 
Eq. JTU). 



Finally, several numerical studies 



41| have focused on the O-point of random bonds models, 
and have argued that its location is extremely well approximated by an annealed computation. 
Once again, this confirms that Eq. (|14|) is a reasonable approximation (the random-bond model is 
recovered by setting |.A| = L, v a = 1/L and E ar i.i.d.'s random variables). This is also related to 
the numerical finding that the global collapse in protein folding dynamics is essentially ^insensitive 
to the specific structure of the sequence, but only depends on its global composition 
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III. GLASS PHASES 

If we follow the entropy density s{(5) of the liquid solution as a function of temperature, we 
find that in any heterogeneous sequence s{j3) turns negative at sufficiently low temperatures. This 
indicates the existence of a phase transition to a glass phase which breaks the translational invari- 
ance. 

As we will show, this glass transition can be of two types. In certain sequences the "entropy 
crisis" is preceded by a local instability of the cavity recursions ©-© around the liquid fixed 
point p* a . This implies the divergence of a properly defined spin-glass susceptibility and signals a 
continuous glass transition towards a phase with fully broken replica symmetry. 

In other sequences, and in the Gaussian random bond model, this local instability is irrelevant 
since it occurs - if at all - in the region of negative entropy of the liquid globule. The glass transition 
is thus necessarily discontinuous (1RSB), as was predicted from replica calculations for the random 
bond model 

Dealing with the glass phases requires some modifications of the simple Bethe Peierls approxima- 
tion which we have been using so far. In this section we will describe first some general properties of 
the glass phases, and explain the general technical tools that can be used to study glass transitions 
using the cavity method. 

A. Proliferation of pure states 

In a glassy phase, the space of conformations is expected to split up in a multitude of pure 
states that are separated by large free energy barriers. The slowest time scale of the system, 
corresponding to jumps between pure states, increases dramatically. 

In mean field approximation, or on the Bethe lattice, this time scale diverges and ergodicity 
is broken at the "dynamic" phase transition. The system eventually undergoes a "static" phase 
transition (with a non-analyticity in the thermodynamic potentials) at a lower temperature 3^| . 

In a finite-dimensional model the "dynamic" phase transition becomes a crossover where the 
nature of the most important dynamical processes changes. Whether the "static" phase transition 
survives in a given model, or not, is not known in general. We shall not enter this dispute here 
since we have little to say about it. In any case, the mean-field-like Bethe approximation, assuming 
the existence of many pure states, yields some useful insight on the glass phase. 

Within one pure state, the conformational probabilities on a given site are well-defined 
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FIG. 3: Schematic phase diagram of copolymers as a function of inverse temperature (3 and chemical 
potential //. A polymer in equilibrium with the solvent is described by the coexistence line. Beyond the 
O-point, /3 > j3& it is in a collapsed phase with a finite density. Depending on the sequence correlations of 
the copolymer there may be a local instability of the liquid (dash-dotted line) , giving rise to a continuous 
glass transition at Pi (see upper graph). In the absence of a local instability down to a critical temperature 
in the range of « 1.23, a discontinuous glass transition will take place. The thermodynamic (static) phase 
transition at (3 a is preceded by a dynamic glass transition at (3d where the phase space splits up into different 
pure states. In the glass phase, the critical chemical potential depends on whether the dynamically relevant 
threshold states (dashed line) or the states dominating the static equilibrium (solid line) are described. 

However, there is no reason to assume the equality of local fields on different sites. Rather one 
expects that in a given pure state the sites will have different preferences for certain polymer 
conformations. 

To proceed, one has to use a statistical description of local fields. We shall not explain here all 
the details of this description, but just give the main definitions and refer the reader to 39, 40] for 
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detailed discussions. In a glassy phase, the number of pure states A/V(w) increases exponentially 
with the volume of the system. The complexity is the monotonously increasing, concave 

function defined by ~ exp(V£(u;)). The natural order parameter is the distribution of local 

fields over the pure states 7 whose free energy density u; 7 is fixed to a value ujq: 

pip) oc J2$(P ~ P (l ' 7) )^K " "*>)■ (15) 

7 

An alternative description consists in using a Legendre transformation of the complexity, by intro- 
ducing the parameter m = (1 / (3)T,' (luq) and working at fixed m instead of fixed ujq 42 1. This com- 
putation is equivalent to a 1RSB scheme with Parisi parameter m. From the free energy at fixed m, 
4>i(m), the complexity is obtained through the Legendre transform: mf3(pi(m) = m(3u) — X(u;). 

In a system with a discontinuous (1RSB) glass transition, this approach gives a full description. 
The complexity is strictly positive in the interval u s < uo < lo^, corresponding to the interval 
mj < m < m s in the 1RSB parameter. The thermodynamically dominant metastable states are 
obtained by minimizing the one-replica free energy lo — In an intermediate temperature 

regime T s < T < T^, the minimum is attained for some free energy cj* (corresponding to m, = 1), 
with lo s < uj* < tod- Below the glass transition, T < T s , the minimum is attained at the lower edge 
uj* = u> s (with S(cj*) = 0), corresponding to the 1RSB parameter < m* < 1. 

In a system with a continuous glass transition (FRSB), the full solution should involve grouping 
states into clusters, and clusters into superclusters, building up a continuous ultrametric hierarchy. 
The approach above amounts to a 1RSB approximation of this full structure, and we shall not 
attempt to go beyond this level of approximation. 



B. Order parameters 

In this section we present two types of order parameters which can be used to identify the glass 
phase. 

For a polymer in Euclidean space, described by the position R{ of monomer i, let us consider 
two replicas of the polymer in the same pure state. In the glass phase, provided the global rotation 
symmetry is broken, the local conformation of the two polymers will have a certain tendency to 
be the same while the liquid phase is completely disordered in this respect. In order to measure 
this effect, we introduce the scalar product of the distance vectors between nearby monomers in 
the replicas (1) and (2): 

% 
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We shall be interested in computing the average of this quantity when the replicas are constrained 
to remain in the same pure state. More precisely, we want to evaluate 

7 

where we average over all states 7 with their Boltzmann weigth w~. This quantity is accessible 
numerically. We consider a polymer which is thermalized at time in a configuration Ri(t = 0). 
We let it evolve for a time t, to a configuration Ri{t). The order parameter is given by the quantity 

(if ,2) ) =1^— [ tMAX dt I - Ri(t)) ■ (R\ + d(o) - Ri(0))) , (is) 

\ /state \tMAx7o N ^ / { ^ =q)} 

evaluated over timescales £max which are large but much smaller than the typical timescale for 
interstate transitions or even full equilibration (in particular much smaller than the time scale for 
diffusion or rotation of the polymer, which diverges with N). 

A simpler order parameter can be defined by first introducing, on each site i of the lattice, the 
quantity Sj which takes the value Si = 1 if the site is occupied by a monomer A, Sj = — 1 if there is 
a B monomer and Sj = if the site is empty. Then the overlap between two configurations 1 and 
2 of the polymer can be defined as 

V 



,(1,2) = 

V 



i=l 

Again, one can compute the typical distance \Qab/ between two conformations in the same 

\ / state 

state by recurring to dynamical simulations. 

fi 2) (1 2) 

Notice that both q\' B and F\ ' define a notion of distance (or similarity) between polymer 

(1 2) 

configurations. However, they describe two complementary aspects of the polymer: q\' B essentially 

characterizes the bias of single sites towards a specific monomer type, whereas the order parameters 
(1 2) 

' measure the conformational similarity of the replicas in the vicinity of a given site, once the 
monomer on that site has been fixed. They measure the freezing of the local degrees of freedom of 
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68l | . In contrast the parameter q^' B 



the polymer's backbone, similarly to the approach of 
is hardly sensitive to the geometric constraints induced by the backbone. 

A dynamical evaluation of the above order parameters is particularly convenient on finite- 
dimensional lattices. Notice that the equilibrium probability for two independent replicas to have 
(1 2) 

a finite overlap q\ U , vanishes with the volume of the lattice because of translation invariance. 

On the Bethe lattice it is more natural to work at a finite monomer density, (see Sec. IV D(l . 
In this case, the random structure of the lattice acts as a "pinning field" , and two replicas of the 
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same system typically have a finite overlap. Following the practice from spin-glass theory, we shall 
measure the probability distribution of the quantity l|T$j) with respect to the Gibbs measure: 

P A B(q) = (5(q-q^ ] )) . (20) 
In a liquid phase, (q^B / vanishes and the function PAs(q) is a ^-function. In a glass 



state 

phase (q^B / > an d the function P^siq) becomes non-trivial, with support in the interval 

\ / state \ / state \ / state \ / state 

vanish at the transition point, while they exhibit a jump in the discontinuous case. 



IV. METHODS TO STUDY THE GLASS PHASES IN THE CAVITY APPROACH 

In this section we present the methods that we use to study the glass transition on the Bethe 
lattice. They are applied to various types of polymers in the next sections. 



A. Local instability towards a soft glass phase 

The simplest glass transition is the one associated with an instability of the liquid. The liquid 
solution is always embedded in the 1RSB formalism as the single pure state that exists at high 
temperature: it is described by the field distribution p(p) = 6(p — p*). This solution becomes 
locally unstable if fluctuations around p* grow on average under the cavity recursion This 
phenomenon occurs when 

k>L* > 1 (21) 

where A max is the largest eigenvalue of the transfer matrix for the propagation of deviations from 
the liquid under the recursion (@J)-(|2|, 

M aa , = dI a [pW, . . . ,p«]/^«| pW=p , . (22) 

(Notice that the stronger instability fc|A max | = 1 13] is irrelevant on a random lattice, since it is 
associated to the establishment of a crystalline order that is inherently frustrated because of the 
presence of large loops.) Beyond the local instability, the distribution of local fields p(p) becomes 
non-trivial, but it remains centered around the unstable liquid fixed point. In physical terms this 
indicates that phase space begins to divide up into a small number of states that comprise a large 
number of microconfigurations. These states are characterized by weak local preferences for certain 
polymer conformations that deviate only slightly from the homogeneous liquid state. 
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The instability (|2Tj) generally develops below a temperature Tj. Calling T cr i s the temperature 
where the entropy vanishes, one can have two types of situations: 

• When Ti < T cr i s , the local instability of the liquid is clearly irrelevant, and a discontinuous 
glass transition must take place at some temperature > T cr i s . 

• When T cr i s < Ti, either the instability drives a continuous glass transition (as we will see 
in specific examples, this seems to be the generic case when the instability occurs in a 
region where the liquid entropy is still large), or there exists again a discontinuous glass 
transition taking place at temperatures T > T and rendering the instability irrelevant. It 
is also possible, that a first continuous glass transition towards a slightly frustrated phase 
undergoes a successive discontinuous phase transition at lower temperatures where a stronger 
degree of freezing takes place. 

Because of the relative simplicity of the liquid phase, it turns out that the stability condition 
(|21j) can be studied explicitly for AB copolymers with an interaction matrix which is symmetric 
under the A <-> B interchange. The detailed calculation is given in Appendix El The dangerous 
eigenvalues A of the matrix M. in (|2*2*)) are found to obey the equation 

± 1 wsinhjP) = A(l - (k\)~ L ) 

kl + wcosh((3) ( k -2) + k(k\)- L + 2(k-l)J2f= l 1 qi(k\)- i ' 

where the sign corresponds to ferromagnetic (+) and antiferromagnetic (-) interactions, respec- 
tively. The temperature dependent parameter w = Yla=iP2a/Po characterizes the liquid solution 
and is independent of L, cf. App. [X]and Eqs. (|B1|) . I)B2|) . The sequence properties enter the above 
expression only through the autocorrelation function qi = (l/L) Yl^=\ a a<Ja+i- 

The local instability /3j occurs at the smallest value of (5 where the characteristic polynomial 
()23|) has a root with \X\ 2 k = 1. Usually, for attractive interactions between equal monomers, 
the relevant eigenvalue is A = \j\fk~ while the instability occurs in general with A = —\j\fk in 
ampholytes. 

The location of the instability for the various types of interactions and sequences will be stud- 
ied in the next sections. Let us just mention here that the (periodic) Gaussian random bond 
heteropolymer generically undergoes a discontinuous 1RSB glass transition, in agreement with 
previous studies [(3(3]. 
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B. Cavity recursion within the 1RSB approximation 

In order to study the glass phase itself, we need to compute the distribution of local fields of 
(|15|) for the Bethe lattice. We shall do it here within the 1RSB cavity formalism of ( jsH l^flj]). 
We shall not rederive the full formalism but give the main ingredients needed for our study. The 
statistical average of the simple cavity recursion (|4I7|) . which holds within a given pure state, leads 
to a recursion relation for this distribution: 

P(P) = ^1 n^P W ) *(P - J [P (1) > • • • 'P^) e-™^ (1) <-<P (fc) ] (24) 
J i=l 

where /[p^ 1 ^, . . . , p^ 8 '] is given by (|l )l -l(7 |l . and Z is a normalization. The non trivial reweighting, 
which depends on the parameter m defined in Section Till Al involves the free energy change induced 
by the recursion, which is given by A/[{pW}] = -/9 -1 log(C7[{pW}]) , where C[{p«}] is the 
normalization term appearing in (|4" |) -l|7 j) . This reweighting accounts for the fact that the number 
of pure states increases exponentially with their free energy. 

The free energy is obtained by properly weighting the contributions of different pure states: 

yndp( P w)<({p«}) +^iog |n^(p w x(p (i) 'P (: 



0m(f)i(m) = —log 



(25) 



i=l 

where w s and wi are the site and link partition functions defined in Eqs. and IJ12J) . The 
complexity S(w) is obtained from (j)i{m) through a Legendre transform: m(3(j)i{m) = mfiuj — E(cj). 
Note that the recursion relation ()24l) is the saddle point equation for the functional 4>i(m) with 
respect to p(p). 

Close to a continuous glass transition, p is strongly peaked around the liquid fixed point p*, 
and we can expand the free energy as a function of the moments of the fluctuations p — p* over the 
pure states, as outlined in Appendix El To leading order the corrections to the liquid free energy 
arise from fluctuations in the "replicon" mode, the unstable direction of the transfer matrix (|22|) . 
whose magnitude grows as (Tj — T) 1 / 2 . The continuous glass transition is found to be of third 
order, 



k. ~ A* = c j2^nT ^ ~ Tf + ° ( (T * " ^ (26) 

where c is a positive constant. This is in contrast to discontinuous glass transitions which are 
(generally) of second order in the free energy. 
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V. TWO EXEMPLARY CASES: THE ALTERNATING AMPHOLYTE AND HP MODEL 

In this section we apply the cavity 1RSB formalism to two specific sequences: the regularly al- 
ternating copolymer chains ABABAB . . . for ampholytic and symmetrized-HP interactions. These 
turn out to be rather extreme representatives in the ensemble of possible neutral copolymers, but 
they are the simplest ones, and they exhibit the characteristics of the continuous (ampholyte) and 
discontinuous (HP) transition in a very clear manner. 

The folding of an alternating copolymer on a regular Bethe lattice is a frustrated problem, 
while, clearly, on a regular cubic lattice, it would just behave as a homopolymer with homogeneous 
interactions Eab = e. However, we expect that as soon as a certain number of defects are intro- 
duced in such sequences, their folding on the cubic lattice will be similarly frustrated. In terms of 
Markovian sequences, we consider here the case of ir -C 1. 

While these sequences are expected to behave differently from the alternating one ir = on 
the cubic lattice, it is reasonable to assume that the ir — > limit is smooth on the Bethe lattice. 
Then the Bethe approximation of it <C 1 sequences can be studied using the perfectly alternating 
sequence, as we do here here. Alternating chains are more easily studied with the cavity method, 
since the number of local fields may be reduced to 5 (with 4 independent degrees of freedom): due 
to the inversion symmetry the local conformations reduce to a £ {0, 1A, IB, 2A, 2B}, where 1A 
(IB) comprises the two conformations j A and J. A (f B and j B). The cavity recursion relation 
(|24|) can thus be handled relatively easily, using a population dynamics algorithm described in 
App.El 

A. Ordered structures, correlations, frustration and the order of the glass transition 

Before embarking on the details of the cavity computation for the alternating chains, we present 
here some simple arguments explaining the very different physical nature of the glass phase in the 
alternating ampholyte, which has a continuous transition, and in the symmetrized HP model, which 
has a discontinuous transition. 

Instead of a Bethe lattice, let us consider a regular tree and ask for a maximally dense polymer 
configuration such that all interactions are satisfied (AB interactions in ampholytes and AA or BB 
interactions in the symmetrized HP model). In Fig. 0] we show typical configurations for each case. 
While there is a stratified order in ampholytic configurations that manifests itself in strong long 
range correlations, the symmetrized HP model has an "ordered" structure that is highly correlated 
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FIG. 4: Unfrustrated, maximally dense structures on a tree (k = 3). The ampholyte (left) has an evi- 
dent stratified order and long range correlations of site occupancies. The location of A's and -B's in the 
HP-polymer is correlated with the backbone configuration (thick edges) which makes the distribution of 
monomers look random. 

with the backbone configuration. No long range correlations may persist, and this dense ground 
state is difficult to distinguish from a dense liquid configuration. 

If we turn back to a Bethe lattice, frustration is induced by the presence of large loops. Odd 
loops are inherently frustrated in the ampholyte since they necessarily have to break up the long 
range correlations of the layered structure. This is not the case in the HP-like model where 
most constraints from loops can be satisfied when the backbone is arranged in the right way. In 
other words, the information about local conformations and the associated constraints cannot be 
propagated far away in the case of the HP-like chain, since the correlations of ordered structures 
die out quickly with distance. As long as the density is not too large and there are sufficient voids 
in the globule, a global frustration will not be able to establish. For the ampholyte, however, it 
will be favorable, even at lower density, to develop local (site) preferences for a certain monomer 
type and thus increase the probability of satisfied interactions. This mechanism is at the basis of 
the instability of local fields in the liquid. Note that in the first place this instability is related 
to the type of monomer accommodated on a given site rather than the backbone structure. The 
latter will only come into play at larger densities/lower temperatures. 

This qualitative discussion applies equally to correlated sequences which are not perfectly al- 
ternating but have a strong tendency to alternate (small ir). At the other extreme, if one considers 
the case of tt close to one, where consecutive monomers tend to be alike, one can apply the same 
type of considerations, but with the roles of ampholyte and HP-like chain reversed. We can thus 
conclude that the local instability of a HP-like chain with long blocks of like monomers is associ- 
ated to the appearance of pure states characterized by the same monomer preferences for small 
regions on the lattice. This is reminiscent of the microphase separation (MPS) [18| which has been 
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much discussed in this context and becomes relevant for sequences with a distinct block structure 
2^], However one should remember that the present formulation of the cavity method, 



m 
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which neglects small loops in the lattice, does not allow any quantitative study of this phenomenon 
(this could be addressed using more refined cluster variational methods). 

Repeating the above arguments for more general cases of short range correlated sequences, one 
sees that in general a local instability is favored by sequences whose monomer distribution tends to 
be annealed (e.g., ampholytes with a tendency towards charge alternation along the sequence). It 
is interes ting to note that such 'annealed sequences' naturally result from common protein design 



schemes 



itmg to note tna 



B. The continuous transition in the AB ampholyte 



We start our quantitative study with the alternating ampholyte on a lattice with connectivity 
k + 1 = 6. For this polymer, the local instability of the liquid found from Q23|) develops at the 
inverse temperature ~ 0.7947, much smaller than in most other neutral sequences. The Parisi 
parameter m remains small throughout this phase. 

A closer analysis of the instability shows that the most unstable eigenvector is antisymmetric 
with respect to the exchange of A and B. This indicates that the pure states are essentially 
characterized by the preference of the sites to accommodate one of the two monomer species, in 
agreement with our qualitative discussion. 

On lowering the temperature, the preference of sites for certain conformations (and not only 
for the respective monomers), increases. This could be interpreted as a growing degree of freezing 
that affects larger and larger length scales. 

Figure |S] shows the basic thermodynamic observables p, s, u in the glass phase, computed in 
the cavity method and compares them to the values found in the unstable liquid solution. The 
data have been computed on the coexistence line, i.e., by fixing p c (P) such that the glass static 
free energy vanishes, <f>\{m s \ fj, c ) = (as explained in App. \C$). The strong frustration of the 
polymer can clearly be seen from the suppression of the density in the glass phase that saturates 
around p = 0.71, while in a liquid phase it would tend to p = 1. The entropy crisis of the liquid is 
prevented, the internal entropy of the pure states remaining rather large even at low temperature. 
There is no sign of a strong (discontinuous) freezing transition. 

In App. [E]we explain how to compute the order parameter (|17|) within the cavity approximation. 
The result for the alternating ampholyte is shown in Fig. H3 which again shows a continuous 
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FIG. 5: Alternating ampholyte on a lattice with 6 neighbors per site. The thick lines show the density 
p, entropy s and internal energy u computed in the glass phase using the 1RSB approximation. The thin 
lines give the corresponding values in the liquid solution, which is unstable beyond fli « 0.7947. The glass 
transition is continuous. 
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persistent part of the distance vector between monomers at a distance d — 1,2,3,4 in the backbone (see 



QHl ), plotted versus the inverse temperature (3. 
transition. 



C. The discontinuous transition in the alternating HP model 



The case of the symmetrized-HP alternating sequence, always on a lattice with connectivity 
k + 1 = 6, is extreme in the opposite sense. The liquid solution is always locally stable, even in 
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FIG. 7: Alternating HP-like polymer on a lattice with 6 neighbors per site: The thick lines show the density 
p, entropy s and internal energy u computed in the glass phase using the 1RSB approximation. The thin 
lines give the corresponding values in the liquid solution, which is always locally stable. The glass transition 
is a discontinuous one; it is an almost perfect freezing transition as in the REM. 

the region of negative entropy. However, running the population dynamics algorithm for the 1RSB 
cavity method, one finds a discontinuous glass transition. The dynamic transition takes place at 
Pd « 1.387, just before the entropy crisis of the liquid (/3 C ris = 1.4525). The static phase transition 
follows at (3 S ~ 1.442, in a region of very high density, p ~ 0.95, and almost vanishing entropy. In 
Fig. [3 we plot the density, entropy and internal energy for the alternating HP-polymer along the 
coexistence curve. The internal entropy of the statically dominating pure states is seen to nearly 
vanish in the frozen phase, and the system barely evolves upon lowering the temperature. This 
scenario is very similar to the abrupt freezing encountered in the random energy model (REM). 

The computation of the order parameter (|17|) proceeds as in the case of the ampholyte. The 
result is shown in Fig. |H] and shows clearly the discontinuous transition. 



As we already stressed, one advantage of our approach consists in the possibility of checking 
mean field computations using numerical simulations of well defined polymer models on a Bethe 
lattice. Here we want to demonstrate this feature by considering the alternating AB ampholyte. 

We made extensive simulations on Bethe lattices with connectivity (k + 1) = 6 and volumes V 
ranging from 100 to 800. For all of the data presented in this Section, we fixed (3 = 2.0 above the 



D. Numerical simulations 
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FIG. 8: Alternating HP-like polymer: The structural order parameters (-Fj ) for the glass phase, 

\ / state 

plotted versus the inverse temperature (3. Note the considerably higher values than in Fig. indicating a 
much stronger freezing of local conformational degrees of freedom. 
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FIG. 9: Average length of the polymers simulated on the Bethe lattice. Sizes of the lattice are indicated in 
the legend. The arrow signal the liquid-soft glass phase transition. 



0-point inverse temperature f3@ ~ 0.693 and varied the chemical potential [x. As ^ is increased, 
the system undergoes at first a second order collapse transition (at /j, ~ —3.21887) and then a 
continuous glass transition to the soft-glass phase (ni{(3) ~ —2.38431). 

Notice that most of the algorithms for simulating polymers on finite-dimensional graphs cannot 
be applied to the Bethe lattice. In fact local moves are impossible because of the absence of short 
loops. On the other hand, global moves would require a detailed knowledge of the loop structure 
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FIG. 10: Simulations of the alternating ampholyte on the Bethe lattice with connectivity k + 1 = 6. The 
energy per site of the polymer (left) and its density (right) are plotted versus the chemical potential. In the 
main frames numerical data (symbols) for various lattice sizes are compared with the cavity results (dashed 
line) for average polymer length (I) — ► oo. The agreement is very good. In the insets we plot the liquid 
prediction for infinite (continuous line) and finite (dashed line) average polymer length, which shows that 
the finite length corrections are already small. Notice that in the density inset the finite-length theoretical 
curve is barely visible because it is superimposed on the data. The arrows indicate the analytic result for 
the glass transition point jii. 



for any graph realization. 

This problem can be overcome by simulating a melt of variable-length polymers, the length 
being finite in the thermodynamic limit. The single-polymer physics is recovered when the average 
length diverges. We refer to App. [F] for a detailed description of our algorithm. In Fig. we 
show our numerical data for the average polymer length (I). Notice that (I) ~ 10 4- 25 within the 
dense phase. As will be clear from the other numerical results, this is enough for assuring small 
deviations from the infinite-length limit. The main effects are: a rounding of the collapse transition 
and a small shift of the soft glass transition (which occurs at finite I) ~ —2.40923). 

In order to achieve equilibration within the soft glass phase we adopted the parallel tempering 
. We tested equilibration using the method of Ref. .7], and always checked the 
acceptance rate for temperature-exchange moves to be larger than 50%. 

In Fig. EJ] we plot the energy per lattice site and the monomer density, as functions of the 
chemical potential fx. Notice that the liquid - soft glass phase transition is barely discernible from 
the monomer density, and the energy curve is also quite smooth. The 1RSB cavity result gives a 
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FIG. 11: Left: Finite-size scaling of the spin glass susceptibility, that exhibits a clear divergence as a function 
of system size at the expected value [i c . Right: probability distribution of the AB overlap, evaluated in the 
soft glass phase at fi = — 1. In the inset of the left frame we compare the second moment of the AB overlap 
(Qab) w ith the analytical prediction. 

very good quantitative description of the transition. 

In order to get a finer description of the glass phase, we have measured the order parameter 
function Pab(q) defined in (|2Uj) . In Fig. ^^ we report our numerical data for this quantity at the 
highest chemical potential considered (fx = —1). Because of the large finite- ^ effects, it would be 
difficult to conclude from the numerics alone that the infinite- F function is non-trivial. However, 
the data agree with the 1RSB predictions for the Edwards- Anderson parameter, <7ea ~ 0.259. 

In the same figure (left frame) we consider the spin-glass susceptibility: 

XSG = ^E^ = V ( [<1ab( S V,s(V)] 2 ) , (27) 

i,j 

This quantity diverges as fi — > fi~ in the thermodynamic limit. In a finite size sample, its behavior 
is ruled by the usual finite-size scaling form 

XSG(V,^ = V 2 ^x[V 1/u ^-^)]. (28) 

From the cavity solution of the model, one finds that (q\ B ) ~ A{f3)\pL — fi c ((3)] 2 for /j, > fi c {P)- 
This result implies the following relation between the critical exponents denned in Eq. (|28j) : 



2-77 + 2/z/ = 1. 



(29) 
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FIG. 12: Average AB overlap {c[ab) among two replicas coupled through a term of the type —NeqAB- The 
full line is the prediction from the cavity method for the zero e limit. Here V — 10 4 . 

In fact we find a nice collapse of data corresponding to different sizes using v = 4 and rj = 3/2. 
The comparison of (q\ B ) with the 1RSB cavity prediction is quite good. 

An alternative approach for exploring the low energy structure of the system consists in cou- 
pling two replicas through their overlap, cf. Eq. H19j) . In practice, one adds a term of the form 
—N (5eqAB{s^- l \ s^ 1 ) to the two-replica Hamiltonian and tries to estimate c/ea as follows 

q EA = lim Jim <(?ab(s«, s&)) Ne . (30) 

In Fig. ^]we show the numerical results for (^ab) jv e on a l ar g e si ze lattice (V = 10 ) and several 
values of e. In order to simulate large lattices, we did not use parallel tempering here. Furthermore, 
we adopted a weaker equilibration criterium, requiring {qab)n e t° be roughly time-independent 
on a logarithmic scale. Once again, the numerical results compare favorably with the outcome of 
the cavity calculation. 

VI. RANDOM MARKOVIAN COPOLYMERS 

One can show using the formula (|23|) that the local instability appears the earlier, the stronger 
the tendency of monomers to be annealed along the sequence, that is, the more ^4's and -B's tend 
to alternate in an ampholyte, or to form blocks in an HP model. In both cases the autocorrela- 
tion function is large and its sign oscillates (alternating sequence) or remain positive ('blocky' 
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sequence) . 

To be more quantitative, let us consider a random copolymer chain in the limit L — > oo char- 
acterized by the probability ir £ [0, 1] of two neighboring monomers to be of the same type. The 
autocorrelation function of such a chain is (in the L — > oo limit) % = (2tt — 1)*. 

In Figs. EH an dE]we plot the inverse temperature $ at the local instability as a function of 
the parameter tt for the ampholyte and symmetrized-HP models. 

This instability is certainly irrelevant when is larger than the inverse temperature of the 
entropy crisis of the liquid, /3 cr i s = 1.4525. This situation occurs for tt > 0.4480 in ampholytes, and 
proves the existence of a discontinuous transition. But already when is smaller than, but close 
to, /3 cr is) one should expect a discontinuous 1RSB transition to take place at a f3 < fy. 

In order to complete the diagram, we have numerically solved the cavity recursion by population 
dynamics for neutral sequences of period L = 20, but otherwise random composition. From the 
experience gained for the extreme case of the alternating HP-model (see below), we expected a 
kind of frozen solution with rather strong local conformational preferences to dominate the low 
temperature phase. Such a solution is rather non-trivial to find in a huge functional space, in 
particular since it has to be expected that it occurs in a discontinuous manner and cannot in 
general be found by randomly perturbing the liquid solution. 

We therefore proceeded by initializing the population in a highly polarized state that we will 
discuss in more detail in the next Section. This state actually corresponds to an unstable fixed 
point, but it turns out that at low temperatures, it is usually quite close to a stable non-trivial 
solution of the 1RSB cavity equations. At each temperature, we iterated the cavity recursion for 
about 100 sweeps of the population dynamics, cf. App. [01 fixing the chemical potential to its 
liquid critical value, since this value describes correctly the thermodynamic equilibrium up to the 
static phase transition. The Parisi reweighting parameter was set to m = 1 in order to detect 
the dynamic transition, i.e., the local instability of the frozen solution. For reasons of numerical 
stability, we restricted ourselves to sequences with an anti-palindromic structure, i.e., sequences 
invariant under inversion and subsequent exchange of A's and -B's. The field distributions p(p) 
inherit this invariance, and thus in each update of a new cavity field we can decide at random to 
apply a symmetry operation to the new fields first. This stabilizes the iteration since it counteracts 
the numerical tendency to spontaneously break the balance between j- and j-states. Indeed, there 
is a gauge degree of freedom associated to the relative weight of the two orientations of the chain, 
and in general it is difficult to maintain them balanced, while it can be enforced in sequences with 
a palindromic symmetry. The reason to choose anti palindromic rather than palindromic ones is 
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FIG. 13: Phase diagram for ampholytes as a function of sequence correlations and inverse temperature. The 
continuous line indicates the local instability (3i of the liquid as a function of the Markov chain parameter 
it. The points with error bars indicate the dynamic phase transitions (3d found numerically for several 
sequences of period L — 20 whereby we associated an effective parameter 7r c ff to each chain such that 
the local instability predicted from 7r c g- coincides with the actual one. Almost independently of the chain 
composition we find a highly frozen phase beyond (3d ~ 1-23 that is reached via a discontinuous glass 
transition well before the liquid would undergo an entropy crisis at /3 cr i s - For tt < 0.50, this freezing is 
preceded by a continuous glass transition, as predicted from the local stability analysis of the liquid. The 
actual thermodynamic freezing transition occurs at a lower temperature (3 S > (3d- The horizontal lines for 
the static and dynamic transitions are an educated guess for the location of these transitions in the limit 
L — > oo. 

to avoid at the same time an asymmetry between A- and -B-states which likely occurs in small 
populations, in particular in the case of attractive interactions among equal monomers. 

Our findings for the sequences of period L = 20 are summarized in the plots 1131 1141 and 1151 
Figure [Tol shows the variance (square of the standard deviation) of the local field for f (a = 1) over 
the distribution p(p) for several sequences as a function of inverse temperature. This is a measure 
for the degree of the local bias away from the liquid. Almost independently of the particular 
sequence statistics we find that for (3 > (3d ~ 1-23 a strongly frozen phase (with very low internal 
entropy) exists with an associated dynamic transition at (3d- Depending on the sequence statistics, 
the regime of higher temperatures is either entirely liquid (e.g., for it < 0.50 in the ampholytes), 
or exhibits a weaker form of frustration in a phase of presumably fully broken replica symmetry. 
The latter continuously joins the liquid solution at the local instability predicted by (|23|) . For the 
phase diagram in the (3 — //-plane for either of the two scenarios we refer to Figs. |31 
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FIG. 14: The same general phenomenology as in the ampholyte, cf. Fig. holds for the HP-type models, 
but here, the continuous transition takes place at tt > 0.5. Notice that the Tr e e window displayed here is 
larger than in Fig. 
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FIG. 15: The variance (^p^ a=1 ^ — {pi( a =i)) 2 p of a selected local field as a function of inverse temperature 
for a variety of anti-palindromic sequences of period L — 20. In general, there is a very distinct discontinuous 
transition around j3d ~ 1.23, that is preceded by a glassy regime with smaller fluctuations in the local fields 
if the sequence has a tendency for anticorrelation in ampholytes (main frame) or correlation in symmetrized 
HP-like chains (inset). The sequences are characterized by their effective Markov chain parameter 7r c ff as in 
Figs. C21 and [H 

The generic picture of a quench in temperature is thus the following: For ampholyte sequences 
with some tendency to alternation or HP-like-sequences with a preference for block formation, 
there is a continuous glass transition whose location depends strongly on the composition of the 
sequence. The corresponding glass phase is characterized by a relatively weak frustration and a 
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rather small number of states that comprise many microconfigurations with some weak local pref- 
erences for certain conformations. This preliminary glass phase undergoes a further discontinuous 
phase transition at a lower temperature fid ~ 1-23 that is almost independent of the sequence 
structure and might be called the effective freezing transition. For sequences with correlations of 
the opposite kind, the freezing transition is the only phase transition and occurs directly from the 
liquid. It is interesting to note that in numerical simulations of the folding dynamics of neutral 
HP-type copolymers, the dynamical glass transition was also found to be essentially independent 
of the sequence 

It is intriguing that the critical parameter of tt separating the FRSB from the 1RSB freezing 
scenario is very close to tt = 1/2 which corresponds to sequences without correlations. This is 
particularly interesting from the point of view of protein folding. The nature of correlations present 
in the amino acid sequences of natural proteins is still a matter of intensive debate. The analysis 
of Pande et al. [51j argues in favor of a tendency for sequences to be annealed, i.e., to exhibit 
positive correlations in the hydrophilicity and anticorrelations in the charge of amino acids, which 
would suggest a bias towards the FRSB freezing scenario for proteins. However, the studies by 



Irback et al. 
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|29|,|30| rather point towards anticorrelations in the HP-type degrees of freedom 
which would favor a scenario with a direct transition from the liquid to the frozen glass. The 
discrepancies of these studies mainly concern the nature of long range correlations while on the 
level of nearest neighbor correlations, the protein sequences appear to be rather random, having 
7r ~ 1/2 with respect to both charge and hydrophobic/hydrophilic degrees of freedom. It would be 
very interesting to understand whether the folding of natural proteins takes advantage from their 
sequences being very close to the critical border between the two scenarios. On the other hand, 
as mentioned earlier, most protein sequence design schemes tend to result in (partially) annealed 
monomer chains which are therefore likely to exhibit the intermediate soft glass phase. 



VII. THE CLOSE-PACKED LIMIT 



In this section we provide a detailed analysis of the frozen phase in the limit of high density. 
We first show the existence of a special 'REM-like' fully polarized solution of the 1RSB cavity 
equations at temperatures below the liquid's entropy crisis. Then we show that this solution is 
stable in the close-packed limit of high densities. 
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A. A fully polarized solution 

There always exists a 'fully polarized' solution to the cavity equation (|24j) which describes pure 
states consisting of essentially one unique frozen polymer configuration. In each such state, a given 
site only admits one specific local conformation. On averaging over the different pure states, the 
given site will be found in conformation a with frequency w a . The local field distributions then 
take the form 

Ppoi(p) = 5>«(Am)5(p - eW), (31) 

a 

where the fields are defined by ef? = S aa i. This distribution solves the cavity equations 
when the frequencies w a (f3, m) coincide with the local fields of a liquid at the renormalized inverse 
temperature j3' = m(3, i.e., w a ((3,m) = p* a {[3' = mj3). The replicated free energy of this fully 
polarized solution is (f)i(f3,m) = (f>n q (mf3). The internal free energy of the corresponding frozen 
states is related to the liquid quantities via f po \(f3, m) = d(rruj)i(f3, m))/dm = u\\ q {m(3) — fipn q (mf3), 
and the complexity of states is found from S po i(/3, m) = s\\ q {m0). As is evident from the nature 
of the pure states, their internal entropy vanishes. 

Let us for a moment postpone the discussion of the relevance of this solution, and first discuss 
its physical interpretation. At each value of (3 we have to maximize 4>i over < m < 1, under the 
condition X > 0. For temperatures above the liquid's entropy crisis, (3 < /3 cr i s , the maximum is 
attained at m = 1 and we have uj g = ujn q . When (3 > /? cr i s , the static glass transition takes place 
and the free energy freezes to u g = uj\\ q {f3 cr i s ) , the Parisi parameter taking the value m s = (3 CT i S /f3. 
So this solution describes a full freezing of the polymer in some isolated specific configurations, 
taking place at (3 = (3 CT - 1S . Notice that this scenario exactly parallels the one found in the REM. 

Our numerical study of the AB-copolymers in their highly frozen phase (beyond (3d ~ 1.23) 
finds a solution p(p) which is close to the form (|31|) . although small deviations persist, and the po- 
larization is not complete. In the particular case of the alternating chain we numerically confirmed 
that the optimal Parisi parameter is well fitted by m s = T/T s on the coexistence line. 

B. Stability analysis and the limit of maximal density 

Up to this point we have not discussed the range of validity of the polarized solution (|31|) . and in 
particular, its stability. Unfortunately, this is a difficult problem, and we only can provide partial 
answers. 
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The basic idea consists in perturbing the Ansatz ()31j) and checking whether the perturbation 
grows under the cavity iteration (|24j) . A simple perturbation consists in adding to (|31|1 some 
'almost polarized' fields with a small total weight. Namely we take a field distribution of the form 

P(p) = (1 - ae)p pol (p) + e ^ w a p a (p) , (32) 

a 

where p a (p) is concentrated on fields p close to e^ a \ In fact, it is more convenient to think of it 
as a distribution over the 'small' fields p = {p a '}a'=£a- Hereafter, we shall use the notation p a (p) 
instead of p a (p). Finally notice that the p a (p) 7 s need not to be normalized. Normalization is 
enforced by the constant a in Eq. (|32j). 

Plugging the nsatz (|32|) into Eq. (f24"|) we get to linear order in e: 

Pa (p) = k ^ p ( a i- ■ - a k\ao) dp ai (q) 6 (p- I[q;a 2 ...a k ]j . (33) 

ai...a k 

Here we distinguished the distribution on the right-hand side, p a {~) from the one on the left-hand 
side p' a (-). In fact we are interested in the stability of the iteration (|24[) and not just in its fixed 
point. Here P(a\ . . . ak\ao) is the probability of finding conformations a\ . . . a& on the k leaves of 
the branch in Fig. ^ constrained to the root being in conformation «o- This must be computed 
within the solution described by Eq. Q31|) and can explicitly be written in terms of the weights 
w a (f3, m). Finally I[q; a 2 . . . ctk] denotes the 'small' components of the cavity iteration: 
{I[q;a 2 ■ ■■a k ]) a = J Q [q,e( Q2 ) . . .e^} for a / a . 

Instead of continuing in full generality, let us consider the example of an alternating F-model in 
the closed-packed limit with Eaa = Ebb = —Eab = — 1 (remember that in this case we found a 
discontinuous phase transition with a highly polarized low temperature phase, cf. Sec. IV C|) . Eqs. 
((Sl^ reduce to 

fc-l . 

PiAiP) = 52fnS( P 2A,P2B) d PlB (q) 5 (p 1B - e^" 1 " 2 ^) + (34) 

n=0 J 

+9aS(pib,P2b) Jdp 2A (q) 5(p 2A - e'^q 1B ) + 9b5(pib,P2b) jdp 2B {q) 5(p 2A - e^qis) , 

P2a(p) = 2$(pib,P2b) JdpiB(q) 5 (pia - ePqiA - e~^q 2 B^j , (35) 

plus two equations obtained by interchanging A and B. Here we used the shorthand 5(x,y) = 
6(x)S(y) and expanded I[q; a 2 . . . a^] in the delta functions to linear order in Q a for o. ^ ol\. The 
weights {/„} and g A/B are given by 

f _ \ ~ A _-^m(fc-l-2n) _ _ ~ 1 

In ~ (2cosh/3m) fc - 1 V n J ' 9A/B ~ 1 + ' 1 j 



34 




FIG. 16: Instability of the completely frozen solution due to rearrangements of the backbone. 

A little thought shows that, after one iteration of Eqs. (|34|) . (|35j) we can set 

Pia{v) = $(P2A,P2b) Pia^ib(pib) +$(pib,P2b) Pia^2a(P2a) , (37) 
P2a{p) = S(p 1b ,P2b) P2A^ia{pia) ■ (38) 

and that the linearized recursions decouple in the three 'sectors' {1A — > IB, IB — ► 1A}, {1A — > 
2A,2B — > IB}, {IB — ► 2B,2A — ► The first sector corresponds to shifts of the chain and 

turns out to be marginally stable (the function I[q; a<i . . . a^] has to be developed to second order 
in q). 

The other two sectors correspond to structural rearrangements of the backbone and become 
unstable when m/3 < {m(3) c = y c = 1/2 • log(2fc — 3). This instability has a simple physical 
interpretation. The pure states described by {m(5) c have a free energy density f c = 1/2. This 
means that on average, a randomly chosen site has one violated neighboring bond, i.e., one neighbor 
occupied by a monomer of the opposite type. It is thus possible to rearrange the backbone of the 
alternating chain without paying energy by opening the chain at the given site and redirecting it 
in the direction of the violated bond, and propagating the rearrangement through the lattice, see 
Fig.Iini 

For k < 6 the instability appears at a smaller value than the liquid entropy crisis, y c < (3 C ris(p 
oo). Thus, at low temperatures, (3 > y c , the thermodynamically relevant close packed states 
are correctly described by the stable polarized solution with m s = (3 cr - ls / (3. In particular, we 
can immediately deduce the ground state energy of Hamiltonian walks of an alternating HP- 
chain on a fixed connectivity random graph from the value of ^iiq(/9 C ris; P oo): This yields 
0.083686,0.120619,0.172602,0.236348 violated bonds per site for k = 3, 4, 5 and 6, respectively. 

A numerical study of the cavity recursion equations at maximal density actually finds, for k < 6, 
a coexistence of the polarized solution with another solution in some intermediate range [y c ,yt\- 
This is a peculiarity of the infinite \x regime, the numerics at finite but large /i suggesting that 
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the polarized solution is unphysical below y t . However, since yt < /3 C ris> the polarized Ansatz still 
correctly describes the low temperature regime. 

What happens away from the [i — > oo limit? The possibility of voids allows for new terms in the 
sum over conformations, cf. Eq. l|33j) . It turns out that the iterations become unstable in the new 
sectors {0 — >\ a, \ a — > 0} and {0 — *\ a, j a — ► 0}. Physically, this means that the presence of voids 
in the lattice always allows for a rearrangement of the polymer configuration in some (perhaps very 
rare) regions, preventing a complete freezing in a single state. Still, at y > yt a stable fixed point 
close to the polarized solution (|3"Tj) exists. 

Let us finally notice that the stability of the polarized solution can be studied within a larger 
2RSB Ansatz |45j. The results coincide with the simplified treatment presented here. These results 
are further confirmed if one studies the behavior of field distributions in the T — > limit following 
Ref. 441. 



C. Exact enumerations on a cube 



In an attempt to verify the 1RSB or even REM-like nature of heteropolymers, Shakhnovich and 
Gutin have exactly enumerated all conformations of fully compact random 27-mers on a 3 x 3 x 3 



cube and calculated the overlap distribution function P(q) as a function of temperature |26j,|6j. 
They interpreted their results in favor of a REM-like scenario where only a small number of states 
dominated the low temperature regime and P{q) exhibited typical features of a discontinuous glass 
phase. In view of our mean field predictions, one would expect to find a different scenario when 
repeating this analysis for copolymers (with a certain amount of sequence correlations) in their 
soft glass phase. 

We first repeated this enumeration study for random ampholytes, and found a P(q) order 
parameter very similar to the random bond case studied originally 61], in agreement with the 
results of [2^]. However the same analysis done for correlated ampholytic sequences with various 
values of tt did not show any clear dependence on it. This absence of evidence can have two 
origins. On the one hand it might be due to the extreme restrictions that full packing imposes on 
the conformations. We have seen above that the fully dense limit is very subtle since physically 
important degrees of freedom, which are found in a system with voids, are artificially suppressed, 
as has been put forward by many authors 0,0,0. On the other hand it seems that these sizes 
are too small to study the true phase space structure of the glass phase. 
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VIII. DISCUSSION AND CONCLUSION 

The cavity method approaches the lattice heteropolymer problem from a new point of view 
in that it analyzes the conformational degrees of freedom of chains with quenched-in sequences. 
Furthermore, this method allows to study the whole temperature range and describes the — 
collapse and the low temperature physics within the same formalism. In this sense we believe it 
provides an interesting new perspective in the analytic studies of heteropolymer folding. 

With this local approach we have studied the frustration effects on a given site of the lattice. 
We find that the decisive features determining the nature of the low temperature physics are the 
short-range correlations in the monomer sequence. Polymers whose monomer distribution along 
the chain tends to be annealed have a proclivity to undergo a continuous glass transition to a soft 
glass phase before the strong freezing transition takes place. In oppositely correlated sequences the 
freezing occurs directly from the liquid phase. A weakly polarized phase with broken ergodicity 
and a high sensitivity to the specific sequence, has also been observed in the extensive numerical 
analysis of the phase diagram for specific hydrophilic/hydrophobic chains |68(], and the qualitative 
differences found between selected sequences indeed reflect the general tendencies that we predict 
from the cavity analysis of the slightly different but closely related HP-like model. 

The temperature of the dynamic transition at which highly frozen pure states appear is almost 
independent of the correlations in the sequence as we found from the numerical solution of the 
1RSB cavity equations. For the time being we do not have a deeper understanding of this finding, 
which is in accordance with numerical observations in the dynamics of copolymer folding. We 
hope to obtain better analytical insight into this phenomenon from a thorough analysis of the 
stability of the highly polarized low temperature states. This would probably also explain why the 
border between the 2RSB freezing scenario with an intermediate soft glass and the scenario of a 
direct transition liquid- frozen glass is so close to the Markov parameter tt = 1/2, corresponding to 
uncorrelated chains. 

It would be interesting to verify the predictions for Markovian chains experimentally (preferably 
with ampholytes where the pair interactions are rather strong). In fact it is possible to fabricate 
Markovian copolymers from a random polymerization process, where tt can be controlled by chang- 
ing the chemical parameters of the solution. Furthermore, it will be very interesting to review the 
studies of sequence correlations in natural proteins in the light of our findings. 

The results of the cavity method are expected to be exact for polymers on random (Bethe) 
lattices, as is indeed corroborated by numerical simulations. However, on real space lattices the 
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Bethe approximation neglects the correlations arising from small loops. 

It would thus be very important to check the effect of sequence correlations through numerical 
simulations of polymers on a cubic lattice, using our mean field predictions as a guideline. One 
regime in which the small loops of the cubic lattice can yield a behavior which is qualitatively 
different from the present mean field analysis is the case where the polymer has a strong tendency 
to form local crumples, as it happens in block copolymers which undergo a microphase separation. 
In order to study such problems analytically, it would be interesting to improve the Bethe approx- 
imation by considering enlarged cavities that contain not only a single site but a small cluster of 
nearby sites. This actually amounts to a further step in the framework of the cluster variational 
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method. For the homopolymeric case a first step in this direction has been carried out in 

Already on the level of the simplest copolymer model we found a surprisingly rich phase diagram 
as a function of temperature and sequence correlations. But clearly, the cavity method is amenable 
to a number of generalizations that allow to study more sophisticated models of biopolymers, 
including for instance backbone stiffness, orientational degrees of freedom, or additional structural 
constraints such as the saturation of monomer-monomer interactions, which are crucial, e.g., for 
the folding of RNA. 

APPENDIX A: FINDING THE LIQUID SOLUTION 

In this Appendix we show how the translation invariant liquid solution can be found by solving 
a set of \A\ + 2 equations (instead of 3L + 1 equations as it may appear from Eqs. ©-(HI)). First 
of all it is convenient to make a change of variables defining 

L (i) («) («) L (i) 

(<) = V ^ z (i) = ^ z (i) = ^ uj« = V ^exnf-/3K ) (Al) 

a=l Po PO PO 0=1 PO 

It is easy to see that the cavity equations ©-(jZI)), the free energy (|TT?|) and all the others observables, 
can be rewritten in terms of these 2L + \A\ + 1 variables. In using the new variables, when not 
specified, we shall assume that the index a belongs to the enlarged space {0} U A. We will set 
Eo, a = -Eo-,0 = 0. 

The liquid fixed point has the translation invariant form Wa = w a , z$ = z-f , z$ = z^ a . The 
corresponding equations are easily written: 

zu = faft-ite" (i±^y" , ( A2) 

1 + w \ l + w J 

h, = k^Sf± (Vr^r' ■ (A3) 

1+wq \ 1 + wo / 
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,r, T - k{k - l)e^y * T ' a+1 ( 1±^V" 2 . (A4) 

a=l 

It is important to notice that the above equations are invariant under the transformation z^ a — > 
7 • zi a , zi a — > 7 _1 • Z[ a for any positive 7: we shall fix this freedom below. The reader can easily 
check that any physical observable (such as the free energy, the local energy or the local density) is 
also invariant under such a transformation. This happens because, when following the chain along 
its conventional direction, each time we arrive at a site i, we are obliged to leave the site as well. 

The above equations admit of course the trivial coil solution z^ a = z^ a = 0. Moreover, if one has 
z^a = (zi ao = 0) for a particular ao, this implies z\ a = (z[ a = 0) for any a. Therefore, we shall 
hereafter assume that z^ a ,z^ a / for any a. In this case Eqs. (|A2|) - (jA3|) imply the consistency 
condition 

Equations (|A2|) and (|A3|) are easily solved: 

*r. = nTT-(TT^) fe "^t. (A6) 



b=a 



ke^ 


f 1 + W Tb 


1 + wo 


V 1 + ^0 


keP» , 


(\ + W ab 


1 + wo 





k-1 



= IIttttIttt?) H> (AT) 
6=1 

where zt, Z\ are two integration constants. We can exploit the invariance mentioned above in order 
to fix z-\ = zi = z. 

Plugging the expressions (|A6|) , (|A7|) into Eq. (|A4|) , and using Eq. (|A5|) , we get 

w a = (k— l)Lz 2 £ ^- . (A8) 

We are therefore left with a set of |.4| + 2 equations (Eq. (jA5|) plus the |.4| + 1 equations in (|A8|0 
for \A\ + 2 real variables (z and the |.A| + 1 variables w a ). As anticipated these equations depend 
on the sequence just through the frequencies u a , a E A. The reader will easily check that the same 
is true for any physical observable. 

Near the point all w a are small, and (|A8|) shows that to lowest order they satisfy w a ~ 
Wq u T e ~^ EaT ■ By imposing that a non-trivial solution of (|A5j) should exist one immediately 

obtains the equation 1)14(1 for the location of the point. 
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APPENDIX B: NEUTRAL AB-COPOLYMERS: LOCAL STABILITY ANALYSIS 

Here we outline the computation of the local stability condition for an AB copolymer having a 
generic period- L sequence. We shall use, depending on the context, the two notations a a G {^4, B}, 
or a a £ {+, — } for the sequence. 

As already mentioned in Sec. IIV Al we consider the case of an interaction matrix symmetric 
under A *-* B interchange. Without loss of generality, we can restrict ourselves to the cases of the 
AF- and F-models defined in Sec. Ill Al Moreover we shall assume that the sequence is neutral, 
i.e. ua, = vb = 1/2. Under these hypothesis, Eqs. (|A2|) - (|A4[) admit the symmetric solution 
z^ a = Z[ a = z/\fL, wq = w, w a = w cosh P, where z and w are determined by solving the equations 

We want to compute the local stability of the cavity recursions (|4* |) -(|7 |) around the above solution. 
We therefore imagine that the cavity fields for one of the sites l,...k (let us say the site 1) have 
been slightly perturbed and compute the effect of such a perturbation on the site 0. To linear order 
we get: 

Sz® = A 6z\% - B Sw™ + C 8wg> , (B3) 
fa<J = A 5^1, - B fa® + C 6w£ , (B4) 

Sw® = Dj2(6z^+6z^)-ESw^+F £ Sw™ , (B5) 

o=l ae{A,B} 
L 

5w® = G ^(e-P^-VSz^ + e- pE "^Sz^l) - H 5w { q ] +F e'^ Sw® . (B6) 

0=1 re{A,B} 



The constants A—H are all positive, and can be expressed in terms of the solution of Eqs. (|B1|) - (|B2}1 . 
In the following we will just need the combinations below: 

A = l, CG= k - r \ W , F= *_l^^^ (B7) 
k ' k 2 L l + cw' 2k 1 + cw ' y J 

where we used the shorthand c = cosh (5. 

We must now identify the most relevant perturbation, i.e., the largest eigenvalue of the linear 

transformation (|B3|) - (jB6() . It is simple to show that the subspace 

L L \ 

8w = o-,y, Sz u = o;Y, 5z ^ = 6w ° = \ ( B8 ) 

0=1 a=l a£{A,B} J 
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is preserved by the iteration ()B3|) - (|B6|) . It can be shown that the most relevant eigenvector lies 
indeed within this subspace. We restrict to it by defining the variables 

L L 
Sw = 5wa ~ Swb , 5-[b = ^2 Sz^ a cr a _ 6 5 [b = ^ $ z la <?a+b , (B9) 

a=l a=l 

where we used a a E {+, — } for the polymer sequence. Using the new variables we can rewrite the 
iteration ()B3|) - (|B6|) as follows: 

4°a = A tu + l + ^Cq- a 6wW, (BIO) 
C = A5fl +1 + ±C qa 5^\ (Bll) 
5w^ = 2Gs{5\\ ) +5$) + 2Fs5w^ , (B12) 

where we introduced the notation s = sinh/3 (for the F-model) or s = — sinh/? (for the AF-model), 
and the sequence correlation function 

1 L 

Qb = j2^2 ° a °a+b ■ ( B13 ) 

a=l 

Notice that q b = q^. This remark allows us to sum Eqs. QBlOj) and (|B11|) and to introduce the 
Fourier transform (for p = 2irn/L, n S {1, . . . , L — 1}): 

L 

S{p)=Y,(^a + h,a)e- tpa . (B14) 

a=l 

We obtain therefore 

<5 (0) (p) = Ae ip <5 (1) (p) + LCq(p) 5w {1) , (B15) 

SwW = 2GsjY, 6il) (P) eip + 2Fs6wil) - ( B16 ) 
p 

We can now set d(°)( p ) = \SW(p), SwW = XSw^, and solve for A, thus recovering Eq. (|23|) . 

APPENDIX C: COEXISTENCE CONDITION FOR A MANY STATES MOLECULE 

It may be interesting to explicitly treat the case of an isolated molecule in equilibrium with 
the solvent and determine the coexistence condition in the glass phase. The result is not obvious 
since the system can exist in many different states 7 £ {1, . . . AT} with (extensive) grand potential 
{r^i, . . . , fl/v - }- Each one of these states describes a molecule confined to a volume V. 
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Let us suppose that each state can be traced as the volume V of the system is changed. This gives 
us the volume-dependent potentials $1^(V). If the state 7 is to describe a molecule in equilibrium 
with the solvent it should exert no pressure on the walls of the container: 

< C1 » 

We want to compute the typical value of the above quantity for states having a certain free- 
energy density: 7 ~ Vuj. Let us step back for a moment and consider the extensive complexity 
V, fi), where we made explicit the dependence upon the volume V and the chemical potential 



fi. If we assume that states do not 
is changed, it is easy to show that 



)ifurcate and do not die (or come into existence) as the volume 



1 



for almost any state 7: 



t^ + dQ^V + dV,n) = S(n 7 ;V,/i). (C2) 

Using the asymptotic behavior £(fi;V,/i) ~ VS(w,/j), and the general relations from Sec. IIIlXI 
we can establish the coexistence condition either in the (m, /X) or in the (uj, fj,) plane (we always 
assume and the energy parameters to be fixed). From (|U2[) . we immediately obtain the condition 
in the (/i, uj) plane: 

(w,yu) = X(u,/i) . (C3) 

This is suggestive of a balance between an "internal" osmotic pressure, u, and an "interstate" 
pressure (E/c^E). In the (m,fi) plane, the condition assumes a more compact form 4>i(m, fi) = 0. 
If we consider the lowest lying states, their free energy density w s (aO is determined by the vanishing 
of the complexity: E(u; s (/x), //) = 0. Therefore Eq. ()C3|) is satisfied for \x = /x s , with w s (fi s ) = 0. 
This coincides with the condition for a unique pure state. If metastable states are considered, 
Eq. (|C3|) receives a non-vanishing contribution from the complexity: in particular, one obtains 
uj > 0. This is quite striking since we did not assume the system to equilibrate among states of 
a given free-energy (which indeed does not happen on the short time scales that are relevant to 
determine the boundary conditions with the solvent). 

In Fig.El we represent the condition ()C3|) in the (u, fj) plane. Notice that in general metastable 
states (with S > 0) on the coexistence line correspond to lower chemical potential than that of 
thermodynamically relevant states. 

Let us finally consider the coexistence line at thermodynamic equilibrium. Dominant states 
are obtained by minimizing the free energy uj — /3 _1 S(w,^) with respect to uj. The coexistence 
chemical potential /i* is then obtained from Eq. ()C3|) . In a more compact (but formal) way, it is 
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CO 








FIG. 17: A schematic view of the coexistence between a multi-state molecule and the solvent. Each line 
represents the evolution of the internal free energy of a state as the volume is changed /i) = const.). 

The thick line shows the states which are in equilibrium with the solvent. In particular, we signal the 
coexistence chemical potentials for static and dynamic states. 

determined from the condition 



In the main body of the paper we focus on the behavior of the polymer on this line. Generally 
speaking, at high temperature the maximum in Eq. (|C4j) is attained at m = 1. Since <pi(m = 
= 0ii q (/i), in this region the coexistence line is the same as for the liquid phase. At lower 
temperatures the maximum is attained for < m, < 1 and the thermodynamic coexistence line 
lies above the liquid one. We refer to Fig. |3]for a summary of this behavior. 

APPENDIX D: EXPANSION OF MOMENTS AT THE CONTINUOUS GLASS 

TRANSITION 

Here we analyze the solution of the cavity recursion near the continuous transition to first 
non-trivial order in an expansion of its moments. 

Using both sides of the cavity recursion equation on the 1RSB level ()24|) in order to calculate 
the moments of the cavity fields, one obtains a set of coupled non-linear equations for the moments 
of the fields p a over the distribution p(p). It is convenient to change coordinates and define the 
fields = A^(p a — p^) hi such a way as to diagonalize the matrix (|22j) . Hereafter we shall 
denote by /j, = 1 the most instable ('replicon') direction in this matrix, and by A the corresponding 
eigenvalue. 




(C4) 
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A careful analysis allows to establish the scaling of the moments with respect to the small 
parameter k\ 2 — 1 ~ Tj — T close to the instability (|21|). The leading moment is given by the 
second moment of the replicon mode. One finds (A 2 ) ~ (Tj — T) (the brackets () denote the 
average with respect to p), while all other moments of deviations from the liquid fixed point p* 
are at least of second order in Tj — T. The only moments of order (T — T) 2 are the first order 
moments of (A^ t ), the remaining second moments (AiA a ), (A^A^), with /i, v > 1 and the higher 
moments (Af), (AjA^) and (Af). 

We will exploit the knowledge about this scaling in order to expand the 1RSB free energy ()25|) 
in powers of T — T around the liquid solution. 

The site term gives rise to a series 



log 



w 



liq 



V \ 1 / \i=l '¥i I 




fc+1 1 fc+1 -. 

^^ M A« + iE^A»A0-) + i 



fc+i 



i=i 



(Dl) 



where summation over direction indices p,, v, p 
the shorthand notation 



1 dlUe 



w 



s OA 



o 2 



w 



s,fj,up 



1 



w s dA^dA^dA^ 



, 3L is tacitly understood and we used 

(D2) 
(D3) 
(D4) 



(i) 

A* J liq 



J liq 



liq 



Note that we have made use of the fact that w s is a multi-linear function of the fields pW so that 
higher derivatives have to occur with respect to variables on different sites. The average () is with 
respect to the distributions p(p^ ),..., p(p^ k+1 ^) on all sites. The link term has an analogous 
expansion as (TmTl . but the triple sum vanishes since there are only two different field variables. 
To proceed, we note the identity 



Ws (pU, . . . , P ( fc+1 )) = c[ P ( 2 ), . . . , P ( fc+1 )]^(p (1) ,/[p (2) , • • • , P (fc+1) ]), 



from which one immediately deduces 



(D5) 



(D6) 
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(2) 

for all directions ji. Using that dC/dA\ J \n q = 0, as follows from the properties of the subspace 



(D7) 



B8j) which the replicon belongs to, one further finds 



and 



Ws,l = wli 



(D8) 



Let us now discuss the terms that appear to increasing order in T — T in the expansion of the 
free energy. There is no first order term proportional to (A 2 ), because of (|D7|) . The second order 
terms (A^), (A^A^), (Af), and (AfA^) come with products of factors w s /i^ and cancel exactly 
between the site and link contributions, due to (jD6|) . The only remaining term to second order is 



(3m (2) 







\ 2 1 





(k + l)k s2,k + l 2 
~ (w s ,n) H — (win) 



(D9) 



However, using (|D8f) . the coefficient in brackets is seen to be of order k\ 2 — 1 ~ T — T. 

The same mechanism suppresses the a priori third order terms (A 2 ) (A^A^) and (A 2 ) (AfA^ 
by an additional factor k\ 2 - 1 while the terms (A?) (A M ), (Af) (Af ) and (Af) (A?) do not 
appear, again because of (|D7|) . The only surviving third order contributions are the site terms 



6 V 2 



proportional to (A 2 ) 3 , 

ftm^ = -(k + l)k(k-l) 
To first non trivial order we finally have 

0i = <hi q + 4> {2) + ^ {3) + o((k\ 2 - 1) 4 ) 
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Since the form (|25|) of the free energy is variational, we may obtain the leading moment to first 
order in T — T from the condition d(pi/d (A|) = 0, 

1 
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(A 2 ) 3 + 0((A;A 2 -1) 



Plugging this result in Eq. (jUl 1|) we get 

(1-m) fc(Jfe + l) 
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Note that the prefactor of (fcA 2 — 1) in Eq. (|D12j) has to be positive for consistency. A negative 
value indicates that there is no stable solution close to the liquid fixed point and the glass transition 
would be discontinuous. By explicit calculation of this coefficient at the instability point we found 
this to happen only in very atypical sequences with highly non-symmetric interactions. 

Evaluating the coefficients ttf Sj ii and w s ,m requires the knowledge of the replicon eigenvector. 
This can be derived for the case of copolymers with symmetric interaction matrix Eaa = Ebb = 
—Eab, and equally frequent monomer species, va = vb = 1/2, extending the arguments of Add.IbI 
In particular we obtain (using the variables defined in Adds. lAl and iBl) : 



oo 



5w = 0, 5z^ a = ^^-^2 i a a+n {k\) n , 5z ia = ^ o a _ n (k\) n , 5w a = a5w . (D14) 

n=0 n=0 

Things simplify considerably in several important cases: (i) alternating copolymers; (ii) anti- 
palindromic sequences; (iii) Markov sequences in the L — > oo limit. In all this cases the ratio 
wJm/w^! vanishes. The basic reason is that, because of Eq. (|D14|) . w Sj ixi turns out to be an 
odd function of {a a }. In these cases the free energy 4>i(m) takes the simpler form, cf. (|26|) . 



(/, ' + 1 ) fc2 1 m ; (k\ 2 -lf + 0{(kX 2 -l) 4 ). (D15) 



12(fc-l) 2 /3 (2-m) 2 

At the glass transition the maximum of (p\ is attained at m s = 0. The fourth order term will 
shift its position to m s tx k\ 2 - 1 ~ Ti — T, as we have explicitly checked in the alternating 
AB- ampholyte . 



APPENDIX E: COMPUTING THE ORDER PARAMETER IN THE CAVITY METHOD 

We show here how to compute the local structural order parameters (|17|) using the cavity 
method. 

In the spirit of the Bethe-Peierls approximation we treat the self-avoidance of the polymer chain 
just on a local level, forbidding it to leave a site on the edge on which it arrived, but neglecting 
further constraints that arise on a real space lattice. In the following, we call "non reversal random 
walks" (NRRW) this restricted class of walks on the cubic lattice. 

Let us rewrite the distance vector between monomers i and % + d as R^ d - = J^t =1 
with ri 1 ^ = Rj^ n — If the positions along the chain are statistically equivalent, the overlap 

state can ^ e written as 

m M = (e (e-+ E ■ (e-+ E *>)) w 

\l=0 \n=l n=l+l / \n=l n=l+l / / state 
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where we split the sum according to the length I over which the replicas stay together and put 
fn = fn = T n for n < I. Note that once I is fixed the common part of the path and the two 
legs of length d — I can be considered as non reversal random walks, only subject to the constraint 
that the legs leave in different directions at the bifurcation. These random walks have all the same 
weight when averaging over pure states. Hence, in order to evaluate (|E1|) it is sufficient to calculate 
the probability (-P(O) state f° r ^ wo re P ncas m the same state to follow the same path over a distance 
I, from which we obtain 

>i 1,2) ) statc = E( P «)state/(M) (E2) 



where 

f(i;d) 




(l>+X>l 2) )) » CE3) 

\n=l n=l+l / I ifDD W m 



NRRW(l) 

the average being taken over the uniform distribution of two NRRW's after I common links. 

Using that in a NRRW one has (f ni ■ r n2 ) NRRW = l/A;'™ 1- ™ 2 ', and distinguishing the different 
possible conformations at the bifurcation, one easily finds 

I— 1 ^ . ^ J d— I j - d— I d—l . 

J=l ni=l n,2=l ni=l ri2=l 

(The first two terms stem from the self overlap of the common part, the term in the middle is the 
cross term between the common part and a leg that continues straight with respect to the common 
part, and the last term is a negative contribution due to two legs leaving in opposite directions.) 

In the liquid state, (-P(O)iiq i s j us t given by the probability that two NRRW's stay together 
over a distance 

fll,(l = 0) = ^j, (E5) 

^ > °> ~ FTTF < E6 » 

Upon injecting (|E4j) . (|E5|) . and (|E6|) in one may verify that (F d ) n = 0. 

In the glass phase, (-P(O)state ^ s mos t easily evaluated as N(l) \P(l) ) , where iV(Z) is the 
number of rooted NRRW's of length I and P(l) is the probability for two replicas to stay on a 
specific path of length I. 

In the Bethe-Peierls approximation the latter can be computed within an enlarged cavity con- 
taining all sites of the path. The average over the states is done by averaging independently over 
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the local field distributions on all neighboring sites, taking into account proper weighting factors: 

/6 m \ 1 J"IW, dp (P M ) Vl;a Wi ;tat ({p t } te T,) m 

\ u /sta tc m^*p(p®)wi**(.{v>}«x) m ' 

where we have introduced the set of indices 2i labeling the neighbors of the I + 1 sites on the path: 

J, = < =1 {(0,i )} Uj,"i (u^W'.t,,)}) Uti {(Ml)} • (E8) 

T-^a denotes the probability, given the local field configuration, for two replicas to both stay on the 
given path up to site I and to separate afterwards, under the condition to start off at site with 
monomer a, 

y*. w^ jl) w!- h) 

v _ ^n±n l;a+ lja+ (Fq \ 

(i) 

The weights are the Boltzmann factors associated with a polymer starting with monomer 

a on site 0, staying on the path, and leaving it at the site I via neighbor 



k JS>j') k \ l-l /fc-i \ f Jlj) k \ 

^•'=i Wa i=i y /'=i \i=i / y w a ±i i=i J 

the sign ± indicating that monomer indices increase/decrease along the path. Notice that in Eq. 
QE9JI we selected arbitrarily one of the two equivalent directions. In the above formula?, W/ ; t t an d 
Wi ;a are the Boltzmann factors associated with the ensemble of all possible configurations on the 
path, and of the configurations restricted to have a monomer a on site 0, respectively. They are 
conveniently calculated recursively via 



w l;a/tot (W^) = c ( P (^ . . . >P <w) w^wtot ({pWiIp^ = 1 (p (u) > • • • >p (W )) ' 

(Ell) 

where / denotes the cavity iteration functional as defined by (jHEJ), and C is the corresponding 



normalization constant. The initial conditions for ()E11|) are simply 

fc+i fc+i 



Wo , a ( P (^) , . . . , P (^+ 1 )) = eft. y: ^a-iV^T n ^ ( ei2 ) 



and 



»W (p (0,1) ' • • • 'P ( °' fc) ) = (p (0 ' 1} ' • • • 'P ( °' fc) ) • (E13) 



48 



APPENDIX F: MONTE CARLO SIMULATIONS ON THE BETHE LATTICE 

In this Appendix we describe our approach to numerical simulations of heteropolymers on the 
Bethe lattice. In the first part we define a model for finite length polymers. In the second one we 
present our Monte Carlo algorithm. 

1. Finite length polymers 

We consider a modified ensemble for a varying number of finite length random walks. More 
precisely, a configuration is defined by n mutually- avoiding SAW's. The chain i shall contain 
Ni monomers, with N\ + . . . + N n = N. The Hamiltonian (^Q) receives contributions both from 
self-contacts within a single chain and from mutual contacts between different chains. The grand- 
canonical free energy is 



- f3u) L (0,(J; Vend 



) = lirn^ ^Eg log [ ^ cndn E e ^ N E ^ mNi - ] ] • (Fl) 

~ X \">0 N>0 w J 

We introduced the chemical potential ^ cn( j which couples to the number of chain ends in the 
solution (or, equivalently, to the number of polymers). The single-polymer ensemble is recovered 
in the ^ e nd — > — oo limit. 

Extending the cavity formalism to the finite-// en d case is quite straightforward. As an illustra- 
tion, we can easily write down the generalization of Eqs. l@J)-(j7J): 

vf = C-'ll4\ (F2) 
i=l 

p (0) = c _ le „ ^ 4j) i^ nd + (F3) 

j=i I i=i J 

p (0) = c-^e^Yl^ L^+j^pfLi) , (F4) 
i=l I i=l J 

p£? = ^V*]^) je 2 ^^ , (F5) 

1=1 ^ 1=1 U7^2 J 

where we used the shorthands pf^-i — Pj.a-i/Va ' ^T»+ 1 ~ P^l+i/^^ ■ 



2. The Monte Carlo algorithm 



As already mentioned in Sec. IV Dl numerical simulations of long fixed-length polymers are 
quite difficult on the Bethe lattice. We thus resort to simulating the variable- length ensemble 
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FIG. 18: The three moves used in our Monte Carlo simulations on the Bethe lattice. The monomers (chain 
links) which change because of the move are represented with hatched circles (wiggly lines). 

corresponding to the free-energy (|Flj) . The algorithm includes three types of moves illustrated 
graphically in Fig. 1181 (a) monomer insertion/deletion; (b) chain extension/reduction; (c) two 
chain junction/disjunction. It is straightforward to show that these three moves ensure ergodicity. 

At each step of the algorithm the type of move and the location in the graph are chosen 
randomly. The move is then accepted according to the Metropolis rule in such a way as to satisfy 
detailed balance with respect to the variable length ensemble (|F1|) . Evidently the algorithm is 
more efficient for moderate lengths of the polymers, i.e., not too large values of |/i e nd|- It can be 
therefore convenient, for producing equilibrated configurations, to gradually decrease |// e nd| to the 
desired value. 

APPENDIX G: NUMERICAL SOLUTION OF THE 1RSB CAVITY EQUATIONS WITH 

POPULATION DYNAMICS 

The cavity recursion equation in the form (|24|) suggests a numerical solution by an iterative 
population dynamics 39]: The distribution of local fields p(p) is represented by a (finite) population 
of fields. An iteration step in the dynamics consists in choosing at random k "parent" members 
pW of the population and calculating the corresponding cavity field p(°) = /({p^}) from @- 
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(J7J). This new field is then exchanged against an old field in the population with probability 
C[{p {i) }] m /C^ x , proportional to the reweighting C[{pW}] m (normalized so as to make sure that 
the probability never exceeds 1). If the dynamics converges to a stationary distribution, its density 
satisfies the recursion equation (|24|l. 

In the soft glass phase, the iteration converges rapidly since the distribution of fields remains 
centered around the unstable liquid fixed point. However, the algorithm considerably slows down 
in the frozen phase where the fields have strong biases towards given conformations. Since the 
biases of the k parent members are only rarely compatible with each other, the reweighting is 
usually very small. The population dynamics is then dominated by rare events with a low degree 
of frustration. Obviously, the probability of frustrated events rapidly increases with the number 
of different local conformations and thus with the length of the period L. For this reason we have 
limited our numerical simulations in the frozen glass phase to populations of 4000 fields for chains 
with L = 20. 
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